module working_valuer_mod
USE functions_mod
USE parameters_mod
IMPLICIT NONE
CONTAINS


subroutine working_valuer(age,capital,temp,tr,r,&
        psi,wage,&
        shock_pi,vprime_temp,v_out,k_out,l_out,partprod,hbar)
INTEGER, intent(in)   :: temp,age
DOUBLE PRECISION,    intent(in)   :: capital,tr,r,psi,partprod,hbar
DOUBLE PRECISION,    intent(in)   :: vprime_temp(shocksN,temp),wage,shock_pi(shocksN)
DOUBLE PRECISION,    intent(out)  :: v_out,k_out,l_out
DOUBLE PRECISION                  :: vprime(temp),l_guess,k_out_temp,v_out_temp,l_out_temp
Integer                           :: e_it

vprime=0.0D0

do e_it =1,shocksN
    vprime=vprime+vprime_temp(e_it,:)*shock_pi(e_it)
end do

    v_out=-brent(0.0D0, 0.925D0, 0.95D0,obj_1,.000000001D0,l_out,capital,age,tr,r,&
        wage,psi,vprime,partprod,hbar,temp,k_out)

l_guess=max(l_out-.02D0,l_out/2.0D0+.01D0,.01D0)
    v_out_temp=-brent(l_guess/2.0D0, l_guess, min(.95D0,l_guess*1.5),obj_1,.000000001D0,l_out_temp,capital,age,tr,r,&
        wage,psi,vprime,partprod,hbar,temp,k_out_temp)
if(v_out_temp>v_out)then
    l_out=l_out_temp
    k_out=k_out_temp
    v_out=v_out_temp
end if

l_guess=min(l_out+.02D0,l_out*1.5D0-.01D0,.95D0)
    v_out_temp=-brent(l_guess/2.0D0, l_guess, min(.95D0,l_guess*1.5),obj_1,.000000001D0,l_out_temp,capital,age,tr,r,&
        wage,psi,vprime,partprod,hbar,temp,k_out_temp)
if(v_out_temp>v_out)then
    l_out=l_out_temp
    k_out=k_out_temp
    v_out=v_out_temp
end if


end subroutine working_valuer




DOUBLE PRECISION FUNCTION brent(ax,bx,cx,func,tol,xmin,capital,age,tr,rate,wage,psi,&
vprime,partprod,hbar,temp,k_out)
USE nrtype; USE nrutil, ONLY : nrerror
IMPLICIT NONE
    INTEGER, intent(in)   :: temp,age
    DOUBLE PRECISION, INTENT(IN)    :: ax,bx,cx,tol,capital,partprod,hbar
    DOUBLE PRECISION, INTENT(IN)    :: tr,rate,psi,vprime(temp),wage
    DOUBLE PRECISION, INTENT(OUT)   :: xmin,k_out
INTEGER, PARAMETER          :: ITMAX=1000
DOUBLE PRECISION, PARAMETER ::CGOLD=0.3819660D0,ZEPS=1.0e-3_sp*epsilon(ax)
INTEGER                     ::iter
DOUBLE PRECISION            :: a,b,d,ebrent,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm
DOUBLE PRECISION :: func
a=min(ax,cx)
b=max(ax,cx)
v=bx
w=v
x=v
ebrent=0.0
fx=func(x,capital,age,tr,rate,wage,psi,&
    vprime,partprod,hbar,temp,k_out)
fv=fx
fw=fx
do iter=1,ITMAX
    xm=0.5D0*(a+b)
    tol1=tol*abs(x)+ZEPS
    tol2=2.0D0*tol1
    if (abs(x-xm) <= (tol2-0.5D0*(b-a))) then
        xmin=x
        brent=fx
        RETURN
    end if
    if (abs(ebrent) > tol1) then
        r=(x-w)*(fx-fv)
        q=(x-v)*(fx-fw)
        p=(x-v)*q-(x-w)*r
        q=2.0D0*(q-r)
        if (q > 0.0D0) p=-p
        q=abs(q)
        etemp=ebrent
        ebrent=d
        if (abs(p) >= abs(0.5D0*q*etemp) .or. &
            p <= q*(a-x) .or. p >= q*(b-x)) then
            ebrent=merge(a-x,b-x, x >= xm )
            d=CGOLD*ebrent
        else
            d=p/q
            u=x+d
            if (u-a < tol2 .or. b-u < tol2) d=sign(tol1,xm-x)
        end if
    else
        ebrent=merge(a-x,b-x, x >= xm )
        d=CGOLD*ebrent
    end if
    u=merge(x+d,x+sign(tol1,d), abs(d) >= tol1 )
    fu=func(u,capital,age,tr,rate,wage,psi,&
        vprime,partprod,hbar,temp,k_out)
    if (fu <= fx) then
        if (u >= x) then
            a=x
        else
            b=x
        end if
        call shft(v,w,x,u)
        call shft(fv,fw,fx,fu)
    else
        if (u < x) then
            a=u
        else
            b=u
        end if
        if (fu <= fw .or. w == x) then
            v=w
            fv=fw
            w=u
            fw=fu
        else if (fu <= fv .or. v == x .or. v == w) then
            v=u
            fv=fu
        end if
    end if
end do
call nrerror('brent: exceed maximum iterations')
CONTAINS
!BL
SUBROUTINE shft(a,b,c,d)
    DOUBLE PRECISION, INTENT(OUT)   :: a
    DOUBLE PRECISION, INTENT(INOUT) :: b,c
    DOUBLE PRECISION, INTENT(IN)    :: d
a=b
b=c
c=d
END SUBROUTINE shft
END FUNCTION brent




DOUBLE PRECISION FUNCTION obj_1(labor,capital,age,tr,rate,wage,psi,&
    vprime,partprod,hbar,temp,k_out)
    INTEGER, intent(in)   :: age,temp
    DOUBLE PRECISION,    intent(in)   :: capital,labor,tr,rate,psi
    DOUBLE PRECISION,    intent(in)   :: vprime(temp),wage,hbar,partprod
    DOUBLE PRECISION, INTENT(OUT)   :: k_out
    INTEGER :: cap_feasible,low_capital
    DOUBLE PRECISION, allocatable, dimension(:) :: c_upper

    allocate(c_upper(temp))
    c_upper=c_upper_fnc2(labor,wage,capital,kgrid(1:temp),temp,tr,rate)
    if(c_upper(1)<ebar+0.001D0) then
        obj_1=1000000.0D0*abs(c_upper(1)-1.002D0)
    else

    cap_feasible=int(MINLOC(c_upper,1,mask=c_upper.gt.ebar+0.001D0))
    low_capital= minloc(vprime,1, mask = vprime.ge.-10000000.0D0)

    obj_1=brent2(kgrid(low_capital),min(kgrid(cap_feasible)*.25D0,kgrid(cap_feasible)*1.25D0), kgrid(cap_feasible),obj_utils,.000000001D0,k_out,capital,labor,&
        age,tr,rate,wage,psi,vprime,partprod,hbar,temp)
    end if
    deallocate(c_upper)
end function obj_1


DOUBLE PRECISION FUNCTION brent2(ax,bx,cx,func,tol,xmin,capital,labor,age,tr,rate,wage,&
    psi,vprime,partprod,hbar,temp)
USE nrtype; USE nrutil, ONLY : nrerror
IMPLICIT NONE
    INTEGER, intent(in)   :: temp,age
    DOUBLE PRECISION, INTENT(IN)    :: ax,bx,cx,tol,capital,partprod,hbar
    DOUBLE PRECISION, INTENT(IN)    :: tr,rate,psi,vprime(temp),wage,labor
    DOUBLE PRECISION, INTENT(OUT)   :: xmin
INTEGER, PARAMETER          :: ITMAX=500
DOUBLE PRECISION, PARAMETER ::CGOLD=0.3819660D0,ZEPS=1.0e-3_sp*epsilon(ax)
INTEGER                     ::iter
DOUBLE PRECISION            :: a,b,d,ebrent2,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm
DOUBLE PRECISION :: func
a=min(ax,cx)
b=max(ax,cx)
v=bx
w=v
x=v
ebrent2=0.0
fx=func(x,capital,labor,age,tr,rate,wage,psi,&
    vprime,partprod,hbar,temp)
fv=fx
fw=fx
do iter=1,ITMAX
    xm=0.5D0*(a+b)
    tol1=tol*abs(x)+ZEPS
    tol2=2.0D0*tol1
    if (abs(x-xm) <= (tol2-0.5D0*(b-a))) then
        xmin=x
        brent2=fx
        RETURN
    end if
    if (abs(ebrent2) > tol1) then
        r=(x-w)*(fx-fv)
        q=(x-v)*(fx-fw)
        p=(x-v)*q-(x-w)*r
        q=2.0D0*(q-r)
        if (q > 0.0D0) p=-p
        q=abs(q)
        etemp=ebrent2
        ebrent2=d
        if (abs(p) >= abs(0.5D0*q*etemp) .or. &
            p <= q*(a-x) .or. p >= q*(b-x)) then
            ebrent2=merge(a-x,b-x, x >= xm )
            d=CGOLD*ebrent2
        else
            d=p/q
            u=x+d
            if (u-a < tol2 .or. b-u < tol2) d=sign(tol1,xm-x)
        end if
    else
        ebrent2=merge(a-x,b-x, x >= xm )
        d=CGOLD*ebrent2
    end if
    u=merge(x+d,x+sign(tol1,d), abs(d) >= tol1 )
    fu=func(u,capital,labor,age,tr,rate,wage,psi,&
        vprime,partprod,hbar,temp)
    if (fu <= fx) then
        if (u >= x) then
            a=x
        else
            b=x
        end if
        call shft(v,w,x,u)
        call shft(fv,fw,fx,fu)
    else
        if (u < x) then
            a=u
        else
            b=u
        end if
        if (fu <= fw .or. w == x) then
            v=w
            fv=fw
            w=u
            fw=fu
        else if (fu <= fv .or. v == x .or. v == w) then
            v=u
            fv=fu
        end if
    end if
end do
call nrerror('brent2: exceed maximum iterations')
CONTAINS
!BL
SUBROUTINE shft(a,b,c,d)
    DOUBLE PRECISION, INTENT(OUT)   :: a
    DOUBLE PRECISION, INTENT(INOUT) :: b,c
    DOUBLE PRECISION, INTENT(IN)    :: d
a=b
b=c
c=d
END SUBROUTINE shft
END FUNCTION brent2







DOUBLE PRECISION FUNCTION obj_utils(kprime,capital,labor,age,tr,rate,wage,psi,&
    vprime,partprod,hbar,temp)
INTEGER, intent(in)   :: age,temp
DOUBLE PRECISION,    intent(in)   :: capital,labor,tr,rate,psi,kprime
DOUBLE PRECISION,    intent(in)   :: vprime(temp),wage,hbar,partprod
DOUBLE PRECISION                  :: con,vprime_out(1),prod,kprime_pass(1),con_tot,eng,dummy

    prod=1.0D0
    if (labor<hbar) then
        prod=partprod
    end if
    con_tot=workcon(labor,capital,kprime,tr,rate,wage*prod)

    dummy=energy_sharer(0.00D0,con_tot*.15D0,con_tot,energy_share_equation,.00001D0,con,con_tot,age)
    eng=(con_tot-con)/pe
    kprime_pass(1)=kprime


    CALL spline_linear_val1(temp,kgrid(1:temp),vprime(1:temp),kprime_pass(1),vprime_out(1),1)
    obj_utils=-valuer(con,eng,labor,psi,vprime_out(1),adult_equivalence(age))

END FUNCTION obj_utils

subroutine spline_linear_val1 ( ndata, tdata, ydata, tval, yval, ndata2 )
      INTEGER :: left, right,i
      INTEGER, intent(in)  ::ndata,ndata2
      DOUBLE PRECISION,  intent(in) :: tdata(ndata),ydata(ndata),tval(ndata2)
      DOUBLE PRECISION,  intent(out) :: yval(ndata2)
      DOUBLE PRECISION  ::  ypval
ypval=0.0D0
yval=0.0D0
do i=1,ndata2
    call rvec_bracket1 ( ndata, tdata, tval(i), left, right )
  ypval = ( ydata(right) - ydata(left) ) / ( tdata(right) - tdata(left) )
  yval(i) = ydata(left) +  ( tval(i) - tdata(left) ) * ypval
end do
  return
end SUBROUTINE spline_linear_val1

subroutine rvec_bracket1 ( n, x, xval, left, right )
  implicit none
      INTEGER, intent(in)  ::n
      INTEGER, intent(out)  ::left,right
      DOUBLE PRECISION, intent(in)   ::xval,x(n)
  integer i
  do i = 2, n - 1

    if ( xval < x(i) ) then
      left = i - 1
      right = i
      return
    end if
   end do
  left = n - 1
  right = n
  return
end subroutine rvec_bracket1


END module working_valuer_mod


